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G.  B.  McFadden,  B.  T.  Murray,  and  S.  R.  Coriell 
National  Institute  of  Standards  and  Technology 
Gaithersburg,  MD  20899 

M.  E.  Glicksman  and  M.  E.  Selleck 
Department  of  Materials  Engineering 
Rensselaer  Polytechnic  Institute 
Troy,  NY  12181 

Abstract 

During  crystal  growth  from  the  melt,  a fundamental  problem  is 
to  understand  the  interaction  of  the  crystal-melt  interface  with  fluid 
flow  in  the  melt.  This  problem  combines  the  complexities  of  the 
Navier-Stokes  equations  for  fluid  flow  with  the  nonlinear  behavior 
of  the  free  boundary  representing  the  crystal-melt  interface.  Some 
progress  has  been  made  by  studying  explicit  flows  that  allow  a base 
state  corresponding  to  a one-dimensional  crystal-melt  interface  with 
solute  and/or  temperature  fields  that  depend  only  on  the  distance 
from  the  interface.  This  allows  the  strength  of  the  interaction  be- 
tween the  flow  and  the  interface  to  be  assessed  by  a linear  stabil- 
ity analysis  of  the  simple  base  state.  The  case  of  a Taylor-Couette 
flow  interacting  with  a cylindrical  crystalline  interface  is  currently 
being  investigated  both  experimentally  and  theoretically.  We  con- 
sider the  changes  in  the  linear  stability  of  this  system  produced  by 
density-driven  convection  generated  by  the  interaction  of  the  density 
gradients  with  the  gravitational  and  centripetal  acceleration. 


1 Introduction 

The  study  of  crystal  growth  from  the  hquid  or  melt  phase  provides  a rich 
source  of  free  boundary  problems,  the  solutions  for  which  are  of  great  prac- 
tical interest.  (Some  descriptions  of  common  crystal  growth  techniques  are 
given  by  Hurle  and  Jakeman  in  [1]  and  references  therein.)  When  growing 
crystals  of  doped  semiconductors  or  metallic  alloys,  the  concentration  of  so- 
lute at  the  freezing  interface  is  of  special  concern  [2].  In  most  appHcations 
it  is  desirable  to  produce  crystals  with  homogeneous  distributions  of  solute 
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throughout  the  crystal,  and  great  care  is  taken  in  the  design  of  the  crys- 
tal growth  apparatus  to  attempt  to  control  the  concentration  and  thermal 
fields  near  the  interface.  In  directional  solidification  from  the  melt,  for  ex- 
ample, an  idealized  furnace  that  is  free  from  imperfections  would  produce  a 
planar  crystal-melt  interface  with  one-dimensional  temperature  and  solute 
fields,  allowing  solute  to  be  incorporated  uniformly  in  the  growing  crystal. 
In  reahty,  even  such  a planar  crystal-melt  interface  is  subject  to  various 
instabihties  [3,  4]  which  can  result  in  segregation  of  solute  at  the  interface 
and  produce  inhomogeneous  distributions  of  solute  in  the  crystal.  In  addi- 
tion to  instabilities  associated  with  the  interface  itself  [5],  under  terrestriaJ 
growth  conditions  it  is  often  difficult  to  avoid  the  occurrence  of  fluid  flow 
in  the  melt  due  to  natural  convection  [6].  Such  flows  are  themselves  able  to 
cause  imdesirable  segregation  of  solute  [7]  and  may  result  in  the  production 
of  inferior  quahty  crystal.  Avoiding  natural  convection  is  one  of  the  main 
motivations  for  developing  the  capability  of  crystal  growth  under  the  mi- 
crogravity conditions  available  in  low  earth  orbit,  where  the  driving  force 
for  natural  convection  is  lower  by  orders  of  magnitude.  In  addition,  such  an 
environment  allows  more  precise  fundamental  experiments  on  interface  dy- 
namics to  be  performed  without  the  complicating  effects  of  buoyancy-driven 
convection  [8]. 

The  study  of  the  interaction  of  fluid  flow  with  a crystal-melt  interface 
is  thus  an  area  of  fundamental  importance  in  materials  science,  but  despite 
much  recent  research  [9]  the  understanding  of  such  interactions  is  fragmen- 
tary. The  general  problem  combines  the  complexities  of  the  Navier-Stokes 
equations  for  the  fluid  flow  in  the  melt  with  the  nonlinear  behavior  of  the 
free  boundary  representing  the  crystal-melt  interface.  Some  progress  has 
been  made  by  studying  exphcit  flows  that  allow  a base  state  corresponding 
to  a one-dimensional  crystal-melt  interface  with  solute  and/or  temperature 
fields  that  depend  only  on  the  distance  from  the  interface.  This  allows  the 


2 


strength  of  the  interaction  between  the  flow  and  the  interface  to  be  assessed 
by  a linear  stability  analysis  of  the  simple  base  state. 

For  example,  one  can  examine  changes  in  the  morphological  stability  [5]  of 
the  interfcice  in  the  presence  of  flow  in  the  melt.  Specific  flows  that  have  been 
considered  in  this  way  include  plane  Couette  flow  [10,  11],  thermosolutal 
convection  [7,  12],  plane  stagnation  flow  [13,  14],  rotating  disk  flow  [15], 
and  the  asymptotic  suction  profile  [16].  One  can  also  examine  changes  in 
the  hydrodynamic  stability  of  a given  flow  that  occur  when  a rigid  bounding 
surface  is  replaced  by  a crystal-melt  interface.  Examples  here  include  the 
instabilities  associated  with  Rayleigh-Benard  convection  [17],  thermosolutal 
convection  [6,  18,  19],  plane  Poiseuille  flow  [20],  the  asymptotic  suction 
profile  [16],  thermally-driven  flow  in  an  annulus  [20,  21],  and  Taylor  Couette 
flow  [22,  23]. 

In  previoiis  work  [22,  23]  we  have  described  the  interaction  of  a Taylor- 
Couette  flow  with  a cylindrical  crystal- melt  interface  under  the  assumption 
that  the  effects  of  buoyancy  can  be  neglected.  This  preliminary  work  is  in 
support  of  experimental  studies  being  conducted  with  succinonitrile;  for  this 
material  the  crystal-melt  interface  is  predicted  to  have  a significant  effect 
on  the  conditions  for  marginal  stabihty  of  the  flow.  Since  large  tempera- 
ture differences  in  the  system  are  capable  of  driving  natural  convection  in 
the  melt  [24,  25],  it  is  desirable  to  include  such  effects  in  our  theoretical 
treatment  as  well.  In  this  paper  we  reformulate  the  problem  to  include 
the  effects  of  buoyancy  when  the  axis  of  the  cylinders  is  aligned  with  the 
direction  of  gravity.  This  generates  a more  complicated  flow  field  in  the 
bcise  state,  which  then  is  subject  to  not  only  centrifugal  instabilities  but 
buoyant  instabilities  as  well  [26,  27].  We  also  include  the  effects  of  density- 
driven  convection  produced  by  the  interaction  between  the  radial  density 
gradient  and  the  centripetal  acceleration  of  the  azimuthal  flow  in  the  base 
state  [28].  Linear  stability  resiilts  are  obtained  numerically  for  a typical 
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FIG.  1.  Schematic  diagram  of  the  crystalline  inner  annulus  (labeled  “S”) 
surrounded  by  the  liquid  phase  (labeled  “L”).  In  the  base  state  the  unper- 
turbed crystal- melt  interface  is  cylindrical,  with  h(z,4>,i)  = Ri. 

experimental  configuration.  For  these  conditions  we  find  that  the  effect  of 

the  contribution  from  the  interaction  of  the  radial  density  gradient  with  the 

centripetal  acceleration  is  not  significant.  The  natural  convection  causes  a 

slight  stabilization  of  the  system;  however,  there  is  a two-fold  increase  in 

the  wavelength  of  the  most  dangerous  disturbance. 

2 Taylor- Couette  Flow  with  Buoyancy 

We  consider  Taylor-Couette  flow  [29,  30]  in  the  presence  of  a crystal-melt 
interface  f = h{z^(f>,i)  (overbars  will  denote  dimensional  quantities).  A 
cylindrical  coordinate  system  (r,  2, 0)  is  used.  The  melt  occupies  the  region 
h{z,  (/>,<)<  r < R2,  and  the  crystal  occupies  the  region  Rq  < r < h{z^  (f>,  t) 
(see  Fig.  1);  in  the  unperturbed  base  state  the  interface  is  an  infinite  cylin- 
der with  h = Ri.  We  also  consider  the  convection-diffusion  equation  for 
heat  transport  as  well.  We  consider  steady  rotation  of  the  system;  the 
outer  cylinder  is  stationary  and  the  inner  cylinder  and  crystal  rotate  with 
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angular  velocity  Cli.  The  temperature  dependence  of  the  density  is  taken 
into  account  in  a Boussinesq  approximation  [31],  which  allows  convective 
effects  due  to  the  interaction  of  the  density  gradient  with  the  gravitational 
and  centripetal  acceleration  terms.  We  £issume  the  axis  of  the  cylinders  is 
aligned  with  the  gravity  vector,  —gz,  where  z is  the  iinit  vector  along  the 
z-axis. 


The  nonlinear  dimensional  governing  equations  in  the  melt  are  the  con- 
tinuity equation. 


du  u \ dv  dw  „ ,,  , 

_ -I- _ 4- = 0,  (la) 

or  r r ocp  oz 


the  momentum  equations. 


Du  dp 
Di  df 


u 


2 dv 


Po—  + S = /<  I ^ I + 


f2  d<f> 


(lb) 


Dv  1 dp 


_ Dw  dp  _ 


(lc) 

(ld) 


and  the  energy  equation 


DT 

— = /cV^T, 

Dt 


where 


D d _ d V d _ d 


(le) 


and 

-2__^  1_5  1 d^  ^ 

df2  f Qf  f2  Qfjj;!  ^^2 

Here  u,  v,  and  w are  the  velocity  components  in  the  f,  <f>,  and  z direc- 
tions, respectively,  p is  the  pressure,  T is  the  temperature  in  the  melt,  p 
is  the  viscosity  coefficient,  k is  the  thermal  diffusivity  in  the  melt,  and 
p = po{l  — a[T  — Te])  is  the  density,  where  po  and  Te  are  reference  densities 
and  temperatures,  and  a is  the  coefficient  of  thermal  expansion.  In  the 
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Boussinesq  approximation  employed  here  [28,  31],  the  temperature  varia- 
tion in  the  density  is  neglected  in  the  advection  terms  on  the  left  hand  sides 
of  the  momentum  equations. 

The  sohd  is  rotating  with  an  azimuthal  velocity  rfii,  and  the  tempera- 
ture in  the  sohd,  Ts,  obeys 

^ + = (If) 

where  ks  is  the  thermal  diffusivity  in  the  crystal. 

At  the  outer  boimdary  r = R2,  have  u = w = v = 0,  and  T = T2, 
where  T2  is  the  (constant)  temperature  imposed  at  the  outer  cy Under;  T2 
is  assumed  to  exceed  the  melting  point  Tm  of  the  material.  At  the  inner 
boundary  r = we  have  Ts  = Tbi  where  Tb  is  the  (constant)  temperature 
imposed  at  the  inner  cyhnder;  Tb  is  assumed  to  he  below  Tb,.  At  the  interface 
f = h{z,  (f>,t)  the  boundary  conditions  are  [6] 


u = u?  = 0, 


(2a) 


V = A(z,  0,  F)Qi, 

f = fs  = fm-  frr^TK, 

_r 

^ di  ^ \ d<f>  d4>  dz  dz 

, (dfs  Idhdfs  dhdfs\ 

isl— -TT-— 


(2b) 

(2c) 

(2d) 


df  h'^  d<p  d(f> 

where  F is  a capillaiy  length,  Tv  is  the  latent  heat  of  fusion  per  unit  volume 
of  crystal,  K is  the  mean  curvature  of  the  interface,  and  and  ks  are 
the  thermal  conductivities  in  the  hquid  and  sohd,  respectively.  We  have 
assumed  equal  densities  of  crystal  and  melt,  and  equal  heat  capacities  in 
each  phase. 
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3 Base  State 


The  equations  admit  an  annular  base  state  with  h{z,<l),t)  = Ri,  u = 0,  v = 
F(f),  w = W{f)^  p = P(f,  f),  T = 0(r),  and  Ts  = 05(r).  Dimensionless 
variables  may  be  introduced  as  follows.  The  length  scale  is  chosen  to  be  the 
melt  gap  width  L = R2  — Ri,  the  time  scale  is  chosen  to  be  where 

1/  = fi/pQ  is  the  kinematic  viscosity,  the  velocity  scale  is  v/  L,  the  pressure 
scale  is  pou'^ and  the  deviation  of  the  temperature  from  its  value  at  the 
unperturbed  interface  is  measured  in  units  of  the  temperature  difference 
across  the  melt,  AT  = T2  — Tm  + fmT/Ri.  The  melt  then  occupies  the 
region  7//(l  — 77)  < r < 1/(1  — tj),  the  unperturbed  interface  is  located  at 
r = Tj/ {1—7}),  and  the  crystal  occupies  the  region  775/(1  —77)  < r < 77/(1  — 77), 
where  77  = R1/R2  and  tjs  = RqI R2  < V-  (Dimensionless  counterparts  to  the 
dimensional  variables  will  lack  overbars.) 

The  resulting  dimensionless  expressions  for  the  bcLse  state  variables  axe 
as  follows.  The  base  azimuthal  velocity  can  be  written  in  the  form 

, 1 Re 77^ 

(1-772)^  r(l - 772)(1  -77)2’ 

where  Re  = L'^Qi/u  is  the  Reynolds  number.  The  solution  for  the  axial 
velocity  that  is  appropriate  to  model  a closed  system  with  no  net  axial 
volume  flux, 


is  given  by 
W{r)  = 


-G 


16(1-77)2 


/‘«2 

/ rW{r)dr  = 0, 


(3) 


-4(4^ -77^) 


2Nln^ 


In  77 


where  G = gaL^AT/u"^  is  the  Grashof  number,  ^ = (1  — T})r,  axid 

^ (1  - 772)(1  - 3772)  - 477^  In 77 

(1  - 772)2  + (1  - 77^*)  ln77 

The  dimensionless  temperature  fields  are  given  by 

ln(^/77) 


0(r)  = 


ln(l/77)’ 


(4a) 
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and 


e (r) 


(4b) 


The  radial  pressure  gradient  balances  the  centrifugal  force,  and  the  axial 
pressure  gradient  is  constant. 


4 Linearized  equations 


We  next  consider  the  dimensionless  equations  hnearized  about  the  steady 
base  state.  We  Fomier  transform  the  axial  and  azimuthal  coordinates,  which 
introduces  the  axial  wavenumber  a and  the  azimuthal  wavenumber  n,  and 
write 


u(r,2,0,t)  ^ 

[ ° 1 

^ u(r)  ^ 

v(r,2,0,  t) 

y(r) 

v(r) 

w{r,z,(i>,  t) 

W{r) 

w{t) 

p(r,2,<^,f) 

— 

P{r,z) 

+ 

p{r) 

T{r,z,<f>,t) 

0(r) 

0(r) 

Ts{r,z,(f>,t) 

Qs{r) 

0s(r) 

h{z,(f),t)  ^ 

^ ^/(l  - V)  ! 

[ h ) 

exp(<Tf  + in<j>  + iaz), 


where  the  perturbation  amphtudes  (quantities  with  hats)  are  assumed  small. 
The  complex  growth  rate  a = <Tr  + io'i  determines  marginal  stabihty:  the 
flow  is  stable  if  <Tr  < 0 for  all  values  of  n and  a.  If  one  sets  D = d/dr^  the 
linearized  equations  in  the  melt  region  7//(l  — r;)  < r < 1/(1  — Tf)  taJce  the 
form 


. 1 » m ^ . . 

Du  H — u H V + law  = 0, 

r r 


inV . . ^ 

ail  ^ u + iaWu  + Dp  = 


(5a) 

(5b) 


D^u  ^ — Du  — 
r 


a r 


. 2.n.'i  , „(l-e0)l/, 

u —V  + 2 V fc). 
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(5c) 


inV  in 

Gv  H V + iaWv  + uDV  H p — 

r T 


D^v  + —Dv  — 


^ 2m  (1  - eQ)V  . 

V + — —u u. 


inV 

GW  H w 4-  iaW w + uDW  + tap  = 

r 


D^w  H — Dw  — 
r 


2 ^ 
+ — 


+ G0, 


inV  A 


(70  H 0 4-  zaW^0  4-  uDQ  = — Z)‘^0  4 — DQ  — 

r Pr  \ r 

and,  in  the  region  ps/{l  — rj)  < r < rj { {I  — p),  one  obtains 


a + ^ 


(5d) 


0 , (5e) 


1 


(705  4"  mRe05  = ( D^Qs  4 — DQs  ~ 


Tt 

4-  — 


05  . 


(5f) 


Here  Pr  = vj kl  is  the  Prandtl  number,  P,  = i//«5,  and  e = aAT. 

The  linearized  boundary  conditions  are  u = t;  = i2;  = 0 = Oatr  = 
1/(1  — T/),  05  = 0 at  r = 775/(1—7;),  and,  at  the interfacer  = = rj/ (1  —77), 


^ = 0,  (5g) 

w + hDW  = 0,  (5h) 

V 4-  hDV  = Keh,  (5i) 

Q + hDe  = es  + hDSs  = -7  (^j) 

- irCh  = [Oe  - qUBs)  • (5k) 


where  7 = {TmT)/{LAT),  C = {uLv)/{kLAT),  and  q = ks/ki. 


5 Numerical  Results 

The  linearized  equations  constitute  an  eigenvalue  problem  from  which  the 
complex  growth  rate  g = Gr  icn  can  be  determined  for  each  choice  of  the 
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remaming  parameters.  Curves  of  marginal  stability  may  be  obtained  by  set- 
ting (Tr  = 0 and  computing  instead  a pair  (Re,  cr,)  as  a function  of  wavenum- 
ber, which  can  be  done  (see  [20])  in  the  manner  suggested  by  Keller  [32]. 
The  numerical  procedure  was  checked  carefully  against  previously  pubHshed 
linear  stability  results  for  related  problems  [22,  27,  28,  33]. 

To  illustrate  the  relative  importance  of  the  eifects  of  buoyancy  on  the 
stability  of  the  system,  we  use  the  material  properties  of  succinonitrile  [21] 
for  a typical  container  geometry.  We  take  R2  = 1.60  cm,  Ri  = 1.11  cm,  and 
Ro  = 0.458  cm,  giving  rj  = 0.690  and  rjs  = 0.286.  We  assume  a temperature 
gradient  in  the  melt  of  2.5  K/cm,  which  produces  a temperature  difference 
across  the  melt  gap  of  AT  = 1.225  K,  giving  C = 455,  e = 9.923  • lO"'*, 
G = 170,  and  7 = 4.9  • 10“^;  for  simplicity  we  set  7 = 0.  The  thermal 
properities  of  the  liquid  and  solid  are  similar  for  succinonitrile,  and  we  take 
the  ratio  of  thermal  conductivities  to  be  g = 1,  with  = Pg  = 22.8.  We 
restrict  our  discussion  to  axisymmetric  disturbances  (n  = 0),  which  are 
expected  to  be  the  most  dangerous  modes  in  this  configuration. 

Marginal  stability  curves  (<7^  = 0)  are  compared  in  Fig.  2 for  four  dif- 
ferent cases.  The  Reynolds  number  Re  of  the  flow  is  plotted  versus  the 
axial  wavenumber  a of  the  disturbance.  The  curves  in  this  plot  are  com- 
puted with  e = 0,  so  that  the  interaction  of  the  density  gradient  with  the 
centripetal  acceleration  is  neglected. 

The  top  two  curves  correspond  to  the  classical  Taylor-Couette  problem 
with  rigid  isothermal  surfaces  boxmding  the  hquid,  and  the  bottom  two 
curves  correspond  to  the  two-phase  problem  with  a crystal-melt  interface 
a.t  f = Ri.  We  first  consider  the  two  cases  without  buoyancy  {G  = 0). 
The  rigid-walled  system  haa  a critical  Reynolds  number  Re  = 35  (see  Ta- 
ble I),  corresponding  to  the  onset  of  the  classical  secondary  flow  consisting 
of  toroidal  Taylor-vortex  cells.  The  crystal-melt  interface  destabilizes  the 
system,  giving  a smaller  critical  Reynolds  number  Re  = 9;  the  axial  wave- 
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length  of  the  most  dangerous  disturbance,  A = 27r/a,  is  shifted  to  larger 
values.  This  effect  has  been  described  in  our  previous  work  [22,  23];  the  size 
of  the  destabilization  is  foimd  to  vary  with  the  Prandtl  number  of  the  melt. 
The  results  for  a rigid-walled  system  actually  provide  the  limiting  values 
for  a conduction-dominated  system  with  Pr  -C  0.  The  destabilization  due 
to  a crystal-melt  interface  becomes  pronounced  for  convection-dominated 
systems  with  Pr  ^ 1.  For  G = 0 the  disturbances  are  stationary  in  both 
cases,  with  <7^  = 0. 


FIG.  2.  Marginal  values  of  the  Reynolds  number.  Re,  versus  the  axial 
wavenumber  of  the  perturbation,  a,  for  axisymmetric  disturbances,  com- 
paring the  effects  of  buoyancy  (G  = 170)  for  the  rigid-wailed  system  (top 
two  curves)  and  for  the  system  with  a crystalline  inner  cylinder  (lower  two 
curves). 

For  G = 170  the  marginal  curves  are  shifted,  as  indicated  by  the  dashed 
curves  in  Fig.  2.  The  rigid-walled  system  is  destabilized  slightly  from 
Re  = 35  to  Re  = 32,  with  a small  increase  in  wavelength.  This  destabiliza- 
tion of  the  Taylor-Couette  flow  by  the  buoyancy  forces  is  consistent  with 
the  trends  reported  previously  [27].  The  toroidal  cells  are  no  longer  station- 
ary, but  drift  downwards  with  an  axial  phase  velocity  equal  to  —a, /a.  The 
effect  of  buoyancy  when  a crystal- melt  interface  is  present  is  qmte  different. 
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however.  There  is  a slight  increase  in  the  criticaJ  Reynolds  number  when 
buoyancy  is  present,  and  there  is  a two-fold  increase  in  the  wavelength  of  the 
most  dangerous  disturbance.  The  secondary  cells  have  a slow  upward  drift 
(cTi  < 0),  but  the  phase  velocity  is  smaller  by  three  orders  of  magnitude. 
This  decrease  in  phase  velocity  is  likely  due  to  the  strong  coupling  between 
the  interface  and  the  flow:  the  interface  deforms  and  accommodates  the 
cellular  flow  structure,  allowing  the  instability  to  occur  at  lower  Reynolds 
numbers.  For  the  cellular  structure  to  translate  in  the  axial  direction,  how- 
ever, the  deformed  interface  must  melt  and  freeze  for  the  wave  to  propagate, 
which  retards  the  dynamics  of  the  process. 


Table  I 


Pr 

G 

€ 

Re 

CTi 

a 

0 

0 

0 

35.2520 

0 

3.143 

0 

0 

0.001 

35.3777 

0 

3.143 

0 

170 

0 

32.3174 

4.6937 

2.967 

0 

170 

0.001 

32.2991 

4.6881 

2.967 

22.8 

0 

0 

9.4518 

0 

1.782 

22.8 

0 

0.001 

9.4508 

0 

1.782 

22.8 

170 

0.0 

11.1517 

-2.7096  -10-3 

0.904 

22.8 

170 

0.001 

11.1509 

-2.7101  -lO-^ 

0.904 

Critical  values  of  the  Reynolds  number  Re,  the  time  constant  <t,-,  and  the 
critical  wavenumber  a for  the  rigid- walled  case  {Pr  = 0)  and  the  case  of  a 
crystalline  inner  annulus  of  succinonitrile  (Pr  = 22.8).  Buoyancy  forces  are 
absent  for  G = 0,  and  the  interaction  of  the  centripetal  acceleration  with 
the  radial  density  is  absent  for  e = 0. 

In  Table  I we  give  values  of  the  critical  Reynolds  numbers,  the  time 
constant  <t,,  and  the  critical  wavenumber  a for  the  cases  shown  in  Fig. 
2,  and  we  also  give  results  for  the  same  cases  computed  with  a non-zero 
contribution  from  the  interaction  of  the  centripetal  acceleration  with  the 
radial  density  gradient,  using  the  value  e = 0.001.  The  changes  brought 
about  by  including  density  variation  in  the  centripetal  acceleration  are  seen 
to  be  insignificant  for  these  cases.  The  larger  change  in  the  stability  of  the 
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system  with  a radial  temperature  gradient  is  clearly  due  to  gravitational 
forces. 
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